;********************************************
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
;load library for  wetbulb_stull
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/heat_stress.ncl"
;The *area_conserve_remap *and *area_conserve_remap_Wrap* have been *deprecated*.
;The* ESMF regrid* <http://www.ncl.ucar.edu/Applications/ESMF.shtml> library is the recommended replacement. It uses a superior methodology.
load "$NCARG_ROOT/lib/ncarg/nclscripts/esmf/ESMF_regridding.ncl"
;********************************************
 begin

;set baseline
 dirPath = "~/source_data/"
 dirPath2 = dirPath+"/outdoor/"

 mo  = (/"ACCESS-CM2","AWI-CM-1-1-MR","BCC-CSM2-MR","CMCC-CM2-SR5","EC-Earth3","GFDL-ESM4","IITM-ESM", \
         "KIOST-ESM","MIROC6","MIROC-ES2L","MPI-ESM1-2-HR","MRI-ESM2-0"/); 12 models
 case = (/"ssp585" /)
 ncase = dimsizes(case);
 nmod = dimsizes(mo)

 syr = 1990; including historical data to identify 1C warming
 fyr = 2099;
 csyr = sprinti("%0.4i",syr);
 cfyr = sprinti("%0.4i",fyr);
 nyr = fyr - syr + 1;
 ndays = nyr*365

;read dimension for 1 deg data
  dimPath=dirPath+"/domain/"
  landfrac = "sftlf_fx_360x180_1deg_gn.nc"
  area = "areacella_fx_360x180_1deg_gn.nc"

  fdim  = addfile(dimPath+landfrac, "r")
  fdim2 = addfile(dimPath+area, "r")
  lat := fdim->lat
  lon := fdim->lon
  nlat = dimsizes(lat);
  nlon = dimsizes(lon);
  garea = fdim2->areacella  ;m^2
  lfrac = fdim->sftlf/100; %
  wgt := garea*1e-6*lfrac ; sqrkm
  wgt@_FillValue = lfrac@_FillValue; make sure FillValue is set correctly for wgt
  copy_VarCoords(garea, wgt)
  printVarSummary(wgt)
  print("global total land area: "+sum(wgt)+" sq km")


  ;--*read population distribution data from "Fig.S2.Population_map_water_rural.ncl"
  ;1) population spending >30mins to access water (rescaled to UNICEF,WHO gross values)
  ;2) rural farmer population data made from 1km SMOD and 1km POP data from EU GHSL
  ;as described in Supplementary Method 2 "Population distribution data"
  dirNc =  dirPath
  filNc2  = "population_data_watercollector_farmer.nc"
  pthNc2  = dirNc+filNc2
  ncdf2 = addfile(pthNc2 ,"r")
  POPR_1deg = ncdf2->POPR_1deg; population engaged in farming
  POPR_dens = ncdf2->POPR_dens
  POPW_1deg = ncdf2->POPW_1deg; population engaged in water collection
  POPW_dens = ncdf2->POPW_dens

  ;--mask global land area by water collection data
  water_area = where(POPW_1deg .gt. 0, wgt, wgt@_FillValue);
  print("Total area of water collection in 2020: "+sum(water_area)*100*1e-6 +" Mha")

  ;--mask global land area by rural population data
  rural_area = where(POPR_1deg .gt. 0, wgt, wgt@_FillValue); rural area
  print("Rural land area in 2020: "+sum(rural_area)*100*1e-6 +" Mha")


  ;---read G0days grided data---
  dirNc  = dirPath2       ; where nc file will be saved
  filNc  = "Gday-gridcell-warminglevel-10yravg-ens-median-vcbc-outdoor.nc"
  pthNc  = dirNc+filNc
  ncdf = addfile(pthNc ,"r")        ; open output netCDF file

  Gday_sun_med  = ncdf->Gday_sun_med; [nlev,nlat,nlon]
  Gday_diff_med  = ncdf->Gday_diff_med; [nlev,nlat,nlon]
  Gday_zero_med  = ncdf->Gday_zero_med; [nlev,nlat,nlon]
  Gday0d_sun_med  = ncdf->Gday0d_sun_med
  Gday0d_diff_med  = ncdf->Gday0d_diff_med
  Gday0d_zero_med  = ncdf->Gday0d_zero_med
  Gday0d_sun_low  = ncdf->Gday0d_sun_low
  Gday0d_diff_low  = ncdf->Gday0d_diff_low
  Gday0d_zero_low  = ncdf->Gday0d_zero_low
  Gday0d_sun_high  = ncdf->Gday0d_sun_high
  Gday0d_diff_high  = ncdf->Gday0d_diff_high
  Gday0d_zero_high  = ncdf->Gday0d_zero_high
  TWday35d_med  = ncdf->TWday35d_med
  TWday35d_low  = ncdf->TWday35d_low
  TWday35d_high  = ncdf->TWday35d_high
  levs2 = ncdf->lev ; 0.8-4.5 by 0.1
  levels = fspan(1,4,31)
  levid = get1Dindex(levs2, levels) ; select 1-4
  nlev = dimsizes(levid)

     
   ;---mask by areas of subpopulation (water collection or farming)
   Amask_3d := conform_dims(dimsizes(Gday0d_sun_med(levid,:,:)),water_area,(/1,2/)) ; conform to 3d: [nlev,nlat,mlon]
   G0days_sun_med1 := where(Amask_3d .gt. 0, Gday0d_sun_med(levid,:,:), Gday0d_sun_med@_FillValue)
   G0days_sun_low1 := where(Amask_3d .gt. 0, Gday0d_sun_low(levid,:,:), Gday0d_sun_low@_FillValue)
   G0days_sun_high1 := where(Amask_3d .gt. 0, Gday0d_sun_high(levid,:,:), Gday0d_sun_high@_FillValue)
   G0days_diff_med1 := where(Amask_3d .gt. 0, Gday0d_diff_med(levid,:,:), Gday0d_diff_med@_FillValue)
   G0days_diff_low1 := where(Amask_3d .gt. 0, Gday0d_diff_low(levid,:,:), Gday0d_diff_low@_FillValue)
   G0days_diff_high1 := where(Amask_3d .gt. 0, Gday0d_diff_high(levid,:,:), Gday0d_diff_high@_FillValue)
   G0days_zero_med1 := where(Amask_3d .gt. 0, Gday0d_zero_med(levid,:,:), Gday0d_zero_med@_FillValue)
   G0days_zero_low1 := where(Amask_3d .gt. 0, Gday0d_zero_low(levid,:,:), Gday0d_zero_low@_FillValue)
   G0days_zero_high1 := where(Amask_3d .gt. 0, Gday0d_zero_high(levid,:,:), Gday0d_zero_high@_FillValue)
   

   ;*area-weighted average number of days per year with Gday>0 in impacted grid cells
   ;only consider grid cells with Gday>0 under each scenario and masked by specific populations
   wgt_sun := where(Gday_sun_med(levid,:,:) .gt. 0, Amask_3d, Amask_3d@_FillValue); weights of those grid cells concerned
   G0_avgdays_sun_med1 = tofloat(dim_sum_n(G0days_sun_med1*wgt_sun, (/1,2/))/dim_sum_n(wgt_sun, (/1,2/)))
   G0_avgdays_sun_low1 = tofloat(dim_sum_n(G0days_sun_low1*wgt_sun, (/1,2/))/dim_sum_n(wgt_sun, (/1,2/)))
   G0_avgdays_sun_high1 = tofloat(dim_sum_n(G0days_sun_high1*wgt_sun, (/1,2/))/dim_sum_n(wgt_sun, (/1,2/)))
   wgt_diff := where(Gday_diff_med(levid,:,:) .gt. 0, Amask_3d, Amask_3d@_FillValue); weights of those grid cells concerned
   G0_avgdays_diff_med1 = tofloat(dim_sum_n(G0days_diff_med1*wgt_diff, (/1,2/))/dim_sum_n(wgt_diff, (/1,2/)))
   G0_avgdays_diff_low1 = tofloat(dim_sum_n(G0days_diff_low1*wgt_diff, (/1,2/))/dim_sum_n(wgt_diff, (/1,2/)))
   G0_avgdays_diff_high1 = tofloat(dim_sum_n(G0days_diff_high1*wgt_diff, (/1,2/))/dim_sum_n(wgt_diff, (/1,2/)))
   wgt_zero := where(Gday_zero_med(levid,:,:) .gt. 0, Amask_3d, Amask_3d@_FillValue); weights of those grid cells concerned
   G0_avgdays_zero_med1 = tofloat(dim_sum_n(G0days_zero_med1*wgt_zero, (/1,2/))/dim_sum_n(wgt_zero, (/1,2/)))
   G0_avgdays_zero_low1 = tofloat(dim_sum_n(G0days_zero_low1*wgt_zero, (/1,2/))/dim_sum_n(wgt_zero, (/1,2/)))
   G0_avgdays_zero_high1 = tofloat(dim_sum_n(G0days_zero_high1*wgt_zero, (/1,2/))/dim_sum_n(wgt_zero, (/1,2/)))

   ;--area-weighted average is similar to the simple average bacause most tropical and low latitude grid cells are of similar size
   count_sun1 = new(nlev, "integer")
   count_diff1 = new(nlev, "integer")
   count_zero1 = new(nlev, "integer")
   do ii=0,nlev-1
     count_sun1(ii) = num(.not.ismissing(wgt_sun(ii,:,:))); count number of true values
     count_diff1(ii) = num(.not.ismissing(wgt_diff(ii,:,:)))
     count_zero1(ii) = num(.not.ismissing(wgt_zero(ii,:,:)))
   end do

   ;----mask by rural population
   Amask_3d := conform_dims(dimsizes(Gday0d_sun_med(levid,:,:)),rural_area,(/1,2/)) ; conform to 3d: [nlev,nlat,mlon]
   G0days_sun_med2 := where(Amask_3d .gt. 0, Gday0d_sun_med(levid,:,:), Gday0d_sun_med@_FillValue)
   G0days_sun_low2 := where(Amask_3d .gt. 0, Gday0d_sun_low(levid,:,:), Gday0d_sun_low@_FillValue)
   G0days_sun_high2 := where(Amask_3d .gt. 0, Gday0d_sun_high(levid,:,:), Gday0d_sun_high@_FillValue)
   G0days_diff_med2 := where(Amask_3d .gt. 0, Gday0d_diff_med(levid,:,:), Gday0d_diff_med@_FillValue)
   G0days_diff_low2 := where(Amask_3d .gt. 0, Gday0d_diff_low(levid,:,:), Gday0d_diff_low@_FillValue)
   G0days_diff_high2 := where(Amask_3d .gt. 0, Gday0d_diff_high(levid,:,:), Gday0d_diff_high@_FillValue)
   G0days_zero_med2 := where(Amask_3d .gt. 0, Gday0d_zero_med(levid,:,:), Gday0d_zero_med@_FillValue)
   G0days_zero_low2 := where(Amask_3d .gt. 0, Gday0d_zero_low(levid,:,:), Gday0d_zero_low@_FillValue)
   G0days_zero_high2 := where(Amask_3d .gt. 0, Gday0d_zero_high(levid,:,:), Gday0d_zero_high@_FillValue)

   ;*area-weighted average number of days per year with Gday>0 in impacted grid cells
   ;only consider grid cells with Gday>0 under each scenario and masked by specific populations
   wgt_sun := where(Gday_sun_med(levid,:,:) .gt. 0, Amask_3d, Amask_3d@_FillValue); weights of those grid cells concerned
   G0_avgdays_sun_med2 = tofloat(dim_sum_n(G0days_sun_med2*wgt_sun, (/1,2/))/dim_sum_n(wgt_sun, (/1,2/)))
   G0_avgdays_sun_low2 = tofloat(dim_sum_n(G0days_sun_low2*wgt_sun, (/1,2/))/dim_sum_n(wgt_sun, (/1,2/)))
   G0_avgdays_sun_high2 = tofloat(dim_sum_n(G0days_sun_high2*wgt_sun, (/1,2/))/dim_sum_n(wgt_sun, (/1,2/)))
   wgt_diff := where(Gday_diff_med(levid,:,:) .gt. 0, Amask_3d, Amask_3d@_FillValue); weights of those grid cells concerned
   G0_avgdays_diff_med2 = tofloat(dim_sum_n(G0days_diff_med2*wgt_diff, (/1,2/))/dim_sum_n(wgt_diff, (/1,2/)))
   G0_avgdays_diff_low2 = tofloat(dim_sum_n(G0days_diff_low2*wgt_diff, (/1,2/))/dim_sum_n(wgt_diff, (/1,2/)))
   G0_avgdays_diff_high2 = tofloat(dim_sum_n(G0days_diff_high2*wgt_diff, (/1,2/))/dim_sum_n(wgt_diff, (/1,2/)))
   wgt_zero := where(Gday_zero_med(levid,:,:) .gt. 0, Amask_3d, Amask_3d@_FillValue); weights of those grid cells concerned
   G0_avgdays_zero_med2 = tofloat(dim_sum_n(G0days_zero_med2*wgt_zero, (/1,2/))/dim_sum_n(wgt_zero, (/1,2/)))
   G0_avgdays_zero_low2 = tofloat(dim_sum_n(G0days_zero_low2*wgt_zero, (/1,2/))/dim_sum_n(wgt_zero, (/1,2/)))
   G0_avgdays_zero_high2 = tofloat(dim_sum_n(G0days_zero_high2*wgt_zero, (/1,2/))/dim_sum_n(wgt_zero, (/1,2/)))

   count_sun2 = new(nlev, "integer")
   count_diff2 = new(nlev, "integer")
   count_zero2 = new(nlev, "integer")
   do ii=0,nlev-1
     count_sun2(ii) = num(.not.ismissing(wgt_sun(ii,:,:))); count number of true values
     count_diff2(ii) = num(.not.ismissing(wgt_diff(ii,:,:)))
     count_zero2(ii) = num(.not.ismissing(wgt_zero(ii,:,:)))
   end do

   ;--1) show ensemble median and 25-75th percentile using grid cell average
   ;--A minimum of three grid cells with ensemble median Gday > 0 are used to compute the area-weighted mean annual number of days 
   G0_avgdays_sun_med1 := where(count_sun1 .gt. 3, G0_avgdays_sun_med1, G0_avgdays_sun_med1@_FillValue);  
   G0_avgdays_sun_low1 := where(count_sun1 .gt. 3, G0_avgdays_sun_low1, G0_avgdays_sun_med1@_FillValue); 
   G0_avgdays_sun_high1 := where(count_sun1 .gt. 3, G0_avgdays_sun_high1, G0_avgdays_sun_med1@_FillValue); 
   G0_avgdays_diff_med1 := where(count_diff1 .gt. 3, G0_avgdays_diff_med1, G0_avgdays_diff_med1@_FillValue); 
   G0_avgdays_diff_low1 := where(count_diff1 .gt. 3, G0_avgdays_diff_low1, G0_avgdays_diff_med1@_FillValue); 
   G0_avgdays_diff_high1 := where(count_diff1 .gt. 3, G0_avgdays_diff_high1, G0_avgdays_diff_med1@_FillValue); 
   G0_avgdays_zero_med1 := where(count_zero1 .gt. 3, G0_avgdays_zero_med1, G0_avgdays_zero_med1@_FillValue); 
   G0_avgdays_zero_low1 := where(count_zero1 .gt. 3, G0_avgdays_zero_low1, G0_avgdays_zero_med1@_FillValue); 
   G0_avgdays_zero_high1 := where(count_zero1 .gt. 3, G0_avgdays_zero_high1, G0_avgdays_zero_med1@_FillValue); 

   G0_avgdays_sun_med2 := where(count_sun1 .gt. 3, G0_avgdays_sun_med2, G0_avgdays_sun_med1@_FillValue);  
   G0_avgdays_sun_low2 := where(count_sun1 .gt. 3, G0_avgdays_sun_low2, G0_avgdays_sun_med1@_FillValue); 
   G0_avgdays_sun_high2 := where(count_sun1 .gt. 3, G0_avgdays_sun_high2, G0_avgdays_sun_med1@_FillValue); 
   G0_avgdays_diff_med2 := where(count_diff1 .gt. 3, G0_avgdays_diff_med2, G0_avgdays_diff_med1@_FillValue); 
   G0_avgdays_diff_low2 := where(count_diff1 .gt. 3, G0_avgdays_diff_low2, G0_avgdays_diff_med1@_FillValue); 
   G0_avgdays_diff_high2 := where(count_diff1 .gt. 3, G0_avgdays_diff_high2, G0_avgdays_diff_med1@_FillValue); 
   G0_avgdays_zero_med2 := where(count_zero1 .gt. 3, G0_avgdays_zero_med2, G0_avgdays_zero_med1@_FillValue); 
   G0_avgdays_zero_low2 := where(count_zero1 .gt. 3, G0_avgdays_zero_low2, G0_avgdays_zero_med1@_FillValue); 
   G0_avgdays_zero_high2 := where(count_zero1 .gt. 3, G0_avgdays_zero_high2, G0_avgdays_zero_med1@_FillValue); 


   lev2 = (/1,2,4/)
   id = get1Dindex(levels, lev2)
   print("number of grid cells with Gday(sun)>0 (masked by water population): ")
   print(" at warming level "+levels(id)+" : "+count_sun1(id))
   print("number of grid cells with Gday(diff)>0 (masked by water population): ")
   print(" at warming level "+levels(id)+" : "+count_diff1(id))
   print("number of grid cells with Gday(zero)>0 (masked by water population): ")
   print(" at warming level "+levels(id)+" : "+count_zero1(id))

   print("annual number of days with Gday(sun)>0 among heat impacted grid cells (masked by water population): ")
   print(" at warming level "+levels(id)+" : "+G0_avgdays_sun_med1(id))
   print("annual number of days with Gday(shade)>0 among heat impacted grid cells (masked by water population): ")
   print(" at warming level "+levels(id)+" : "+G0_avgdays_diff_med1(id))
   print("annual number of days with Gday(zero)>0 among heat impacted grid cells (masked by water population): ")
   print(" at warming level "+levels(id)+" : "+G0_avgdays_zero_med1(id))
   print("")

   print("number of grid cells with Gday(sun)>0 (masked by rural population): ")
   print(" at warming level "+levels(id)+" : "+count_sun2(id))
   print("number of grid cells with Gday(diff)>0 (masked by rural population): ")
   print(" at warming level "+levels(id)+" : "+count_diff2(id))
   print("number of grid cells with Gday(zero)>0 (masked by rural population): ")
   print(" at warming level "+levels(id)+" : "+count_zero2(id))

   print("annual number of days with Gday(sun)>0 among heat impacted grid cells (masked by rural population): ")
   print(" at warming level "+levels(id)+" : "+G0_avgdays_sun_med2(id))
   print("annual number of days with Gday(shade)>0 among heat impacted grid cells (masked by rural population): ")
   print(" at warming level "+levels(id)+" : "+G0_avgdays_diff_med2(id))
   print("annual number of days with Gday(zero)>0 among heat impacted grid cells (masked by rural population): ")
   print(" at warming level "+levels(id)+" : "+G0_avgdays_zero_med2(id))
   print("")
   

   print("distribution of G0days_sun_med2 at 4C warming") 
   G0days_sun_med2 := where(Gday_sun_med(levid,:,:) .gt. 0, G0days_sun_med2, G0days_sun_med2@_FillValue); 
   levid4 = get1Dindex(levels,4)
   opt = True
   opt@PrintStat = True
   cstat = stat_dispersion(G0days_sun_med2(levid4,:,:), opt) ;detailed statistics



  ; define a polygon centered the ensemble median 
  ;================================================
   xx = levels ;show data from 1 to 4 degC by 0.1 degC

   xp    = new( (/2*nlev/), float )
   yp1   = new( (/2*nlev/), float )
   yp2   = new( (/2*nlev/), float )
   yp3   = new( (/2*nlev/), float )
   yp4   = new( (/2*nlev/), float )
   yp5   = new( (/2*nlev/), float )
   yp6   = new( (/2*nlev/), float )
   yp11   = new( (/2*nlev/), float )
   yp22   = new( (/2*nlev/), float )
   yp33   = new( (/2*nlev/), float )
   yp44   = new( (/2*nlev/), float )

   xp(0:(nlev-1)) = xx
   xp(nlev:(2*nlev-1)) = xx(::-1) ; reverse


   ;ensemble CI for G0_days

   ;ensemble CI for grid cells masked by water data
   yp1(0:(nlev-1)) = G0_avgdays_sun_high1 ; uppper bound
   yp1(nlev:(2*nlev-1)) = G0_avgdays_sun_low1(::-1) ; lower bound
   yp2(0:(nlev-1)) = G0_avgdays_diff_high1 ; uppper bound
   yp2(nlev:(2*nlev-1)) = G0_avgdays_diff_low1(::-1) ; lower bound
   yp3(0:(nlev-1)) = G0_avgdays_zero_high1 ; uppper bound
   yp3(nlev:(2*nlev-1)) = G0_avgdays_zero_low1(::-1) ; lower bound

   ;ensemble CI for grid cells masked by rural data
   yp4(0:(nlev-1)) = G0_avgdays_sun_high2 ; uppper bound
   yp4(nlev:(2*nlev-1)) = G0_avgdays_sun_low2(::-1) ; lower bound
   yp5(0:(nlev-1)) = G0_avgdays_diff_high2 ; uppper bound
   yp5(nlev:(2*nlev-1)) = G0_avgdays_diff_low2(::-1) ; lower bound
   yp6(0:(nlev-1)) = G0_avgdays_zero_high2 ; uppper bound
   yp6(nlev:(2*nlev-1)) = G0_avgdays_zero_low2(::-1) ; lower bound


   print("OK")

 ;---draw line plots
 imagePath = dirPath+"/../image/" 
 image2 = imagePath+"/Fig.6"
 wks2 = gsn_open_wks("pdf",image2)

 res = True
 res@gsnDraw              = False
 res@gsnFrame             = False

 ;res@vpXF                  = 0.1;0.15
 ;res@vpWidthF              = 0.65
 ;res@vpHeightF             = 0.50
 res@vpWidthF             = 0.32; 0.6 ;plot width
 res@vpHeightF            = 0.3; 0.3   ;plot length

 res@tiXAxisFontHeightF   = 0.014;0.022
 res@tiYAxisFontHeightF   = 0.014;0.022
 res@tmXBLabelFontHeightF = 0.012
 res@tmYLLabelFontHeightF = 0.012
 res@tmYRLabelFontHeightF = 0.012
 res@gsnStringFontHeightF = 0.014
 res@tmXTOn               = False
 ;res@tmXBMinorOn          = False
 ;res@tmBorderThicknessF    = 1.5

 res@pmLegendDisplayMode    = "Never";"Always"
 res@pmLegendSide           = "Top"; reference axis
 res@pmLegendParallelPosF   = .38; left-right adjustment
 res@pmLegendOrthogonalPosF = -.31;
 res@pmLegendWidthF         = 0.4;0.1
 res@pmLegendHeightF        = 0.04;0.4
 res@lgLabelFontHeightF     = .015
 res@lgPerimOn              = False
 res@lgOrientation = "horizontal"
 res@lgLabelAngleF = 300

 cmap = read_colormap_file("srip_reanalysis");("radar")
 cmap = cmap(::-1,:) ; reverse color map
 cmap := cmap((/1,2,3,4,5,6,7,8,10,11,12,13,15,16,18,0/),:); first black, last gray for ensemble median
 res@tiYAxisString   = "G~S~day~N~ or G~B~TW~N~~S~day~N~ (W m~S~-2~N~)"
 ;res@tiXAxisString   = "Warming since pre-industrial (~S~o~N~C)" 
 res@trXMinF          = 1
 res@trXMaxF          = 4


 ;-----------plot a----------
 res@vpYF          = 0.90   ; location of plot1
 res@vpXF          = 0.10    ; location of plot1

 
  ;----plot a----
  plot1 = new(3,graphic)
  res@tmYLOn              = True
  res@tmYROn              = False
  res@tmYRLabelsOn        = False
  res@trYAxisType = "LinearAxis"
  res@gsnStringFontHeightF = 0.014;
  res@gsnLeftStringOrthogonalPosF = 0.03 ;adjust the space/distance of gsnLeftString/gsnCenterString
  res@gsnCenterStringOrthogonalPosF = 0.02
  res@gsnLeftString       = "~F22~a";
  res@gsnCenterString  = "Focus on people collecting water"
  ;res@tiYAxisString    = " Mean annual number of days ~C~      with G~S~day~N~>0 or T~B~W~N~~S~day~N~>35";
  res@tiYAxisString    = "Days per year with G~S~day~N~>0";
  res@tiYAxisJust      = "CenterCenter"
  ;res@tiYAxisOffsetXF  = -0.01
  res@trYMinF          = 0;
  res@trYMaxF          := 20; days
  ;res@tmYLMode         = "Automatic"; "Manual"
  res@tmYLMode            = "Manual"
  res@tmYLTickStartF      = 0
  res@tmYLTickEndF        = 20
  res@tmYLTickSpacingF    = 5
  res@tmYLMinorPerMajor   = 4
  res@tmYLAutoPrecision = False
  res@tmYLPrecision     = 2; decimal place

  res@xyLineColors      = "slategray"; "deepskyblue2"
  res@xyDashPatterns    = 0; 3
  plot1(0) = gsn_csm_xy(wks2, xx, G0_avgdays_zero_med1, res); zero solar
  ;plot1(0) = gsn_csm_xy(wks2, xx, G0days_zero_med_mean1, res); zero solar

  gsres = True
  gsres@gsFillOpacityF    = 0.1
  gsres@gsFillColor     = "slategray4";                ; color chosen
  dummy3 = gsn_add_polygon (wks2,plot1(0),xp,yp3,gsres);


  res@tiXAxisString    = " "
  res@tiYAxisString    = " "
  res@gsnLeftString    = " ";
  res@gsnCenterString  = " ";
  res@tmXBOn              = False
  res@tmYLOn              = False
  res@tmYLLabelsOn        = False
  res@tmXBLabelsOn        = False

  ;=====G under full radiation===========
  res@xyLineColors          = "purple2";
  plot1(1) = gsn_csm_xy(wks2, xx, G0_avgdays_sun_med1, res);
  ;plot1(1) = gsn_csm_xy(wks2, xx, G0days_sun_med_mean1, res); zero solar
  gsres@gsFillColor       = "purple";                ; color chosen
  dummy1 = gsn_add_polygon (wks2,plot1(1),xp,yp1,gsres) ;Gday99g_sun

  ;==============under only diffuse radiation==========
  res@xyLineColors        = "orange2"
  plot1(2) = gsn_csm_xy(wks2, xx, G0_avgdays_diff_med1, res);
  ;plot1(2) = gsn_csm_xy(wks2, xx, G0days_diff_med_mean1, res); zero solar
  gsres@gsFillColor       = "orange"                    ; color chosen
  dummy2 = gsn_add_polygon (wks2,plot1(2),xp,yp2,gsres) ;Gday99g_diff
 

  ;--add some text
  txres                     = True                 ; text mods desired
  txres@txFontHeightF       = 0.012                ; default size is HUGE!
  txres@txJust              = "CenterLeft"         ; puts text on top of bars
  ;gsn_text(wks1,plot1,"Hottest 1% of land grid cells", 1.2, res@trYMaxF*0.95,txres) ;

  ;--add legend for shade
  lbres                    = True
  lbres@vpWidthF           = 0.15          ; labelbar width
  lbres@vpHeightF          = 0.075;0.1          ; labelbar height
  lbres@lbBoxMajorExtentF  = 0.6           ; smaller value for more space between color boxes
  lbres@lbMonoFillPattern  = True
  lbres@lbPerimOn          = False
  lbres@lbBoxLinesOn       = False
  lbres@lbFillOpacityF     = 0.2; 0.1

  ;legend for lines, more general function gsn_legend_ndc
  lgres                    = True
  lgres@vpWidthF           = 0.125                ; width of legend (NDC)
  lgres@vpHeightF          = 0.075;0.1           ; height of legend (NDC)
  lgres@lgLineThicknessF   = 1.5                 ; line thicknesses
  lgres@lgPerimOn          = False               ; no perimeter
  lgres@lgLabelOffsetF     = 0.08
  lgres@lgLabelFontHeightF = 0.07;          ; font height (also scales with vpWidthF)


  labels := (/"","",""/); no labels for the shade bars
  lbres@vpWidthF           = 0.18          ; labelbar width
  lbres@vpHeightF          = 0.08;0.10          ; labelbar height
  ;lgres@lbFillColors       := (/"deepskyblue","orange","purple"/)
  lbres@lbFillColors       := (/"slategray4","orange","purple"/)
  ;gsn_labelbar_ndc(wks2,3,labels,0.58,0.9,lbres) ; draw left labelbar column
  gsn_labelbar_ndc(wks2,3,labels,0.07,0.9,lbres) ; draw left labelbar column
  lgres@lgDashIndexes      := (/0,0,0/)          ; line patterns
  lgres@vpHeightF          = 0.08              ; height of legend (NDC)
  lgres@vpWidthF           = 0.11              ; width of legend (NDC)
  lgres@lgLabelFontHeightF = 0.07;          ; font height (also scales with vpWidthF)
  ;lgres@lgLineColors       := (/"deepskyblue2","orange2","purple2"/)
  lgres@lgLineColors       := (/"slategray","orange2","purple2"/)
  ;labels = (/"T~B~W~N~~S~day~N~ > 35","G~S~day~N~ > 0 (diffuse)","G~S~day~N~ > 0 (full solar)"/)
  ;labels = (/"G~S~day~N~ > 0 (dark)","G~S~day~N~ > 0 (shade)","G~S~max~N~ > 0 (sun)"/)
  labels := (/"Dark","Shade","Sun"/)
  ;gsn_legend_ndc(wks2,3,labels,0.63,0.9,lgres)
  gsn_legend_ndc(wks2,3,labels,0.12,0.9,lgres)

  draw(plot1)


 ;---------plot b---------
 ;use vpXF to place two plots horizontally
  res@vpXF          = 0.55       ; location of plot2

  plot2 = new(3,graphic)
  res@tmYLOn              = True
  res@tmYLLabelsOn        = True 
  res@tmXBOn              = True
  res@tmXBLabelsOn        = True
  res@tmYROn              = False
  res@tmYRLabelsOn        = False
  res@trYAxisType = "LinearAxis"
  res@gsnStringFontHeightF = 0.014;
  res@gsnLeftStringOrthogonalPosF = 0.03 ;adjust the space/distance of gsnLeftString/gsnCenterString
  res@gsnCenterStringOrthogonalPosF = 0.02
  res@gsnLeftString       = "~F22~b";
  res@gsnCenterString  = "Focus on agricultural workers"
 ; res@tiXAxisString   = "Warming since pre-industrial (~S~o~N~C)" 
  ;res@tiYAxisString    = " Mean annual number of days ~C~      with G~S~day~N~>0 or T~B~W~N~~S~day~N~>35";
  res@tiYAxisString    = "Days per year with G~S~day~N~>0";
  res@tiYAxisJust      = "CenterCenter"
  ;res@tiYAxisOffsetXF  = -0.01
  res@trYMinF          = 0;
  res@trYMaxF          := 20; days
  ;res@tmYLMode         = "Automatic"; "Manual"
  res@tmYLMode            = "Manual"
  res@tmYLTickStartF      = 0
  res@tmYLTickEndF        = 20
  res@tmYLTickSpacingF    = 5
  res@tmYLMinorPerMajor   = 4
  res@tmYLAutoPrecision = False
  res@tmYLPrecision     = 2; decimal place

  res@xyLineColors      = "slategray"; "deepskyblue2"
  res@xyDashPatterns    = 0; 3
  plot2(0) = gsn_csm_xy(wks2, xx, G0_avgdays_zero_med2, res); zero solar
  ;plot2(0) = gsn_csm_xy(wks2, xx, G0days_zero_med_mean2, res); zero solar

  gsres = True
  gsres@gsFillOpacityF    = 0.1
  gsres@gsFillColor     = "slategray4";                ; color chosen
  dummy4 = gsn_add_polygon (wks2,plot2(0),xp,yp6,gsres);


  res@tiXAxisString    = " "
  res@tiYAxisString    = " "
  res@gsnLeftString    = " ";
  res@gsnCenterString  = " ";
  res@tmXBOn              = False
  res@tmYLOn              = False
  res@tmYLLabelsOn        = False
  res@tmXBLabelsOn        = False

  ;=====G under full radiation===========
  res@xyLineColors          = "purple2";
  plot2(1) = gsn_csm_xy(wks2, xx, G0_avgdays_sun_med2, res);
  ;plot2(1) = gsn_csm_xy(wks2, xx, G0days_sun_med_mean2, res); zero solar
  gsres@gsFillColor       = "purple";                ; color chosen
  dummy5 = gsn_add_polygon (wks2,plot2(1),xp,yp4,gsres) ;Gday99g_sun

  ;==============under only diffuse radiation==========
  res@xyLineColors        = "orange2"
  plot2(2) = gsn_csm_xy(wks2, xx, G0_avgdays_diff_med2, res);
  ;plot2(2) = gsn_csm_xy(wks2, xx, G0days_diff_med_mean2, res); zero solar
  gsres@gsFillColor       = "orange"                    ; color chosen
  dummy6 = gsn_add_polygon (wks2,plot2(2),xp,yp5,gsres) ;Gday99g_sun

  draw(plot2)

  print("finish plot 6")






end 
